# Replication script for "The Unintended Consequences of UN Sanctions: 
# A qualitative comparative analysis", Contemporary Security Policy (2022),
# by Katharina L. Meissner & Patrick A. Mello
# https://doi.org/10.1080/13523260.2022.2059226

# Corresponding Author: Patrick A. Mello
# Website: https://patrickmello.com
# E-mail: p.a.mello@vu.nl

# Load Packages
library(QCA); library(SetMethods); library(ggplot2)

# Read Data
unint <- read.csv("CSP Meissner Mello raw data.csv", header = TRUE, row.names = 1, sep = ",")

# Calibration
unint$AR <- calibrate(unint$AR_raw, type = "fuzzy", method = "direct", c(5, -1, -6), logistic = TRUE)
unint$AR <- round(unint$AR, 2)

unint$LS <- calibrate(unint$LS_raw, type = "fuzzy", method = "direct", c(1, 2, 2.5), logistic = TRUE)
unint$LS <- round(unint$LS, 2)

unint$P5 <- calibrate(unint$P5_raw, type = "fuzzy", method = "direct", c(0.3, 0.4, 0.7), logistic = TRUE)
unint$P5 <- round(unint$P5, 2)

unint$EI <- calibrate(unint$EI_raw, type = "fuzzy", method = "direct", c(80, 50, 40), logistic = TRUE)
unint$EI <- round(unint$EI, 2)

unint$LD <- calibrate(unint$LD_raw, type = "fuzzy", method = "direct", c(48, 96, 144), logistic = TRUE)
unint$LD <- round(unint$LD, 2)

unint$UC <- calibrate(unint$UC_raw, type = "fuzzy", method = "direct", c(0.5, 0.75, 1.25), logistic = TRUE)
unint$UC <- round(unint$UC, 2)

# Data frame
UC <- unint[c(7:12)]
write.csv(UC, "CSP Meissner Mello calibrated data.csv")
head(UC)

# Necessary Condition Analysis (Table 3)
QCAfit(UC[, 1:5], UC$UC, necessity = TRUE)
QCAfit(UC[, 1:5], 1 - UC$UC, necessity = TRUE)

# Conservative Solution
TTcs <- truthTable(UC, outcome = "UC", conditions = "AR, LS, P5, EI, LD",
                    complete = FALSE, show.cases = TRUE, 
                    sort.by = c("incl", "n"), n.cut = 1, incl.cut = 0.75)
TTcs #(Table 4)

SOLcs <- minimize(TTcs, outcome = "UC",include = "1", row.dom = FALSE,
                    details = TRUE, show.cases = TRUE, use.tilde = TRUE,
                    all.sol = TRUE, method = "CCubes")
SOLcs 
  
# Parsimonious Solution
TTps <- truthTable(UC, outcome = "UC", conditions = "AR, LS, P5, EI, LD",
                    complete = TRUE, show.cases = TRUE, 
                    sort.by = c("incl", "n"), incl.cut = 0.75)
TTps #(Table A2)

SOLps <- minimize(TTps, outcome = "UC", include = "1, ?", row.dom = TRUE,
                    details = TRUE, show.cases = TRUE, use.tilde = TRUE,
                    method = "CCubes")
SOLps #(Table A3)

# Intermediate Solution
SOLis <- minimize(TTps, outcome = "UC", include = "1, ?", row.dom = TRUE,
                  details = TRUE, show.cases = TRUE, use.tilde = TRUE,
                  exclude = c(1, 5, 9, 17), method = "CCubes")
SOLis #(Table 5)

## Appendix  

# Figure A1: Histograms of Fuzzy-Set Membership Scores
par(mfrow = c(3, 2), lwd = 2)

h.UC <- hist(unint$UC_raw, freq = TRUE, col = "white", ylim = c(0, 8), xlim = c(0,4),
             xlab = "Unintended Consequences (raw)", ylab = "Frequency", font = 2, 
             font.lab = 2, main = "", lty = 1, lwd = 2)

h.AR <- hist(unint$AR_raw, freq = TRUE, col = "white", ylim = c(0, 4), xlim = c(-10,10),
            xlab = "Autocratic Regime (raw)", ylab = "Frequency", font = 2, 
            font.lab = 2, main = "", lty = 1, lwd = 2)

h.EI <- hist(unint$EI_raw, freq = TRUE, col = "white", ylim = c(0, 6), xlim = c(0,150),
             xlab = "Economic Isolation (raw)", ylab = "Frequency", font = 2, 
             font.lab = 2, main = "", lty = 1, lwd = 2)

h.LS <- hist(unint$LS_raw, freq = TRUE, col = "white", ylim = c(0, 6), xlim = c(0,4),
             xlab = "Large Scope (raw)", ylab = "Frequency", font = 2, 
             font.lab = 2, main = "", lty = 1, lwd = 2)

h.LD <- hist(unint$LD_raw, freq = TRUE, col = "white", ylim = c(0, 6), xlim = c(0,300),
             xlab = "Long Duration (raw)", ylab = "Frequency", font = 2, 
             font.lab = 2, main = "", lty = 1, lwd = 2)

h.P5 <- hist(unint$P5_raw, freq = TRUE, col = "white", ylim = c(0, 10), xlim = c(0,1),
             xlab = "P5 Involvement (raw)", ylab = "Frequency", font = 2, 
             font.lab = 2, main = "", lty = 1, lwd = 2)

# Figure A2: Raw Data and Calibrated Fuzzy Sets
par(mfrow = c(3, 2), lwd = 2)

p.UC <- plot(unint$UC_raw, unint$UC, xlab = "Unintended Consequences (raw)",
             ylab = "Unintended Consequences (fuzzy)", font = 2, font.lab = 2) +
  abline(v = 1.0, col = "black", lwd = 0.7) +
  abline(h = 0.5, col = "black", lwd = 0.7)

p.AR <- plot(unint$AR_raw, unint$AR, xlab = "Autocratic Regime (raw)",
             ylab = "Autocratic Regime (fuzzy)", font = 2, font.lab = 2) +
  abline(v = -1, col = "black", lwd = 0.7) +
  abline(h = 0.5, col = "black", lwd = 0.7)

p.EI <- plot(unint$EI_raw, unint$EI, xlab = "Economic Isolation (raw)",
             ylab = "Economic Isolation (fuzzy)", font = 2, font.lab = 2) +
  abline(v = 50, col = "black", lwd = 0.7) +
  abline(h = 0.5, col = "black", lwd = 0.7)

p.LS <- plot(unint$LS_raw, unint$LS, xlab = "Large Scope (raw)",
             ylab = "Large Scope (fuzzy)", font = 2, font.lab = 2) +
  abline(v = 2, col = "black", lwd = 0.7) +
  abline(h = 0.5, col = "black", lwd = 0.7)

p.LD <- plot(unint$LD_raw, unint$LD, xlab = "Long Duration (raw)",
             ylab = "Long Duration (fuzzy)", font = 2, font.lab = 2) +
  abline(v = 96, col = "black", lwd = 0.7) +
  abline(h = 0.5, col = "black", lwd = 0.7)

p.P5 <- plot(unint$P5_raw, unint$P5, xlab = "P5 Involvement (raw)",
             ylab = "P5 Involvement (fuzzy)", font = 2, font.lab = 2) +
  abline(v = 0.4, col = "black", lwd = 0.7) +
  abline(h = 0.5, col = "black", lwd = 0.7)

# Analysis non-outcome
# Complex solution
NTTcs <- truthTable(UC, outcome = "~UC", conditions = "AR, LS, P5, EI, LD",
                   complete = FALSE, show.cases = TRUE, 
                   sort.by = c("incl", "n"), n.cut = 1, incl.cut = 0.75)

NSOLcs <- minimize(NTTcs, outcome = "~UC",include = "1", row.dom = FALSE,
                  details = TRUE, show.cases = TRUE, use.tilde = TRUE,
                  all.sol = TRUE, method = "CCubes")
NSOLcs #(Table A4)

# Robustness test for calibration thresholds

# Calibration ranges for each condition (Table A1)
rob.calibrange(raw.data = unint, calib.data = UC,
               test.cond.raw = "UC_raw", test.cond.calib = "UC",
               test.thresholds = c(0.5, 0.75, 1.25), step = 0.01, max.runs = 100,
               outcome  = "UC", conditions = c("AR", "EI", "LS", "LD", "P5"),
               incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

rob.calibrange(raw.data = unint, calib.data = UC,
               test.cond.raw = "AR_raw", test.cond.calib = "AR",
               test.thresholds = c(5, -1, -6), step = 0.1, max.runs = 100,
               outcome  = "UC", conditions = c("AR", "EI", "LS", "LD", "P5"),
               incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

rob.calibrange(raw.data = unint, calib.data = UC,
               test.cond.raw = "EI_raw", test.cond.calib = "EI",
               test.thresholds = c(80, 50, 40), step = 1, max.runs = 50,
               outcome  = "UC", conditions = c("AR", "EI", "LS", "LD", "P5"),
               incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

rob.calibrange(raw.data = unint, calib.data = UC,
               test.cond.raw = "LS_raw", test.cond.calib = "LS",
               test.thresholds = c(1, 2, 2.5), step = 0.1, max.runs = 20,
               outcome  = "UC", conditions = c("AR", "EI", "LS", "LD", "P5"),
               incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

rob.calibrange(raw.data = unint, calib.data = UC,
               test.cond.raw = "LD_raw", test.cond.calib = "LD",
               test.thresholds = c(48, 96, 144), step = 1, max.runs = 50,
               outcome  = "UC", conditions = c("AR", "EI", "LS", "LD", "P5"),
               incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

rob.calibrange(raw.data = unint, calib.data = UC,
               test.cond.raw = "P5_raw", test.cond.calib = "P5",
               test.thresholds = c(0.3, 0.4, 0.7), step = 0.05, max.runs = 20,
               outcome  = "UC", conditions = c("AR", "EI", "LS", "LD", "P5"),
               incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

# Consistency threshold range and frequency cutoff range
rob.inclrange(data = UC, step = 0.005, max.runs = 100, outcome = "UC",
              conditions = c("AR", "EI", "LS", "LD", "P5"), 
              incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

rob.ncutrange(data = UC, step = 1, max.runs = 10, outcome = "UC",
              conditions = c("AR", "EI", "LS", "LD", "P5"), 
              incl.cut = 0.75, n.cut = 1, include = "?", exclude = c(1, 5, 9, 17))

### Robust core tests (RC)

# Initial solution (IS)
IS <- minimize(TTps, outcome = "UC", include = "1, ?", row.dom = TRUE,
                   details = TRUE, show.cases = TRUE, use.tilde = TRUE,
                   exclude = c(1, 5, 9, 17), method = "CCubes")

# Alternative test solutions (TS)
# TS1: Higher consistency
TTA1 <- truthTable(UC, outcome = "UC", conditions = "AR, LS, P5, EI, LD",
                   complete = TRUE, show.cases = TRUE, 
                   sort.by = c("incl", "n"), incl.cut = 0.80)

TS1 <- minimize(TTA1, outcome = "UC", include = "1, ?", row.dom = TRUE,
               details = TRUE, show.cases = TRUE, use.tilde = TRUE,
               exclude = c(1, 5, 9, 17), method = "CCubes")
TS1

# TS2: Changed calibration, outcome
UC2 <- UC
UC2$UC <- calibrate(unint$UC_raw, type = "fuzzy", method = "direct", c(0, 1, 2), logistic = TRUE)
UC2$UC <- round(UC2$UC, 2)

TTA2 <- truthTable(UC2, outcome = "UC", conditions = "AR, LS, P5, EI, LD",
                   complete = TRUE, show.cases = TRUE, 
                   sort.by = c("incl", "n"), incl.cut = 0.75)
TTA2

TS2 <- minimize(TTA2, outcome = "UC", include = "1, ?", row.dom = TRUE,
                details = TRUE, show.cases = TRUE, use.tilde = TRUE,
                exclude = c(1, 5, 9, 17), method = "CCubes")
TS2

# TS3: Changed calibration, long duration 
UC3 <- UC
UC3$LD <- calibrate(unint$LD_raw, type = "fuzzy", method = "direct", c(60, 120, 180), logistic = TRUE)

TTA3 <- truthTable(UC3, outcome = "UC", conditions = "AR, LS, P5, EI, LD",
                   complete = TRUE, show.cases = TRUE, 
                   sort.by = c("incl", "n"), incl.cut = 0.75)
TTA3

TS3 <- minimize(TTA3, outcome = "UC", include = "1, ?", row.dom = TRUE,
                details = TRUE, show.cases = TRUE, use.tilde = TRUE,
                exclude = c(1, 5, 9, 17), method = "CCubes")

TS3

# Creating test set TS
TS <- list(TS1, TS2, TS3)
rob.corefit(test_sol = TS, initial_sol = IS, outcome = "UC")
rob.fit(test_sol = TS, initial_sol = IS, outcome = "UC")

# Case-oriented robustness
rob.xyplot(test_sol = TS, initial_sol = IS, outcome = "UC",
           all_labels = TRUE, fontsize = 3, jitter = TRUE, area_lab = TRUE)

RCCtest <- rob.cases(test_sol = TS, initial_sol = IS, outcome = "UC")
RCCtest$CaseParameters
RCCtest$CaseTypes
rob.singletest(test_sol = TS, initial_sol = IS, outcome = "UC")

### Correlations between indicators of the outcome (Table A5)
outcome <- read.csv("CSP Meissner Mello correlations.csv")
outd <- outcome[, 2:6]

options(digits=3)
cor(outcome[,2:6])

library(corrgram)
corrgram(outcome, order = TRUE, lower.panel = panel.shade, upper.panel = panel.pie)

library(psych)
corr.test <- corr.test(outcome[, 2:6], use = "complete")
corr.test

library("Hmisc")
outd.rcorr <- rcorr(as.matrix(outd))
outd.rcorr

outd.coeff <- outd.rcorr$r

library(corrplot)
corrplot(outd.coeff)

### End